%***********************************************************************************
% Code that generate results in the paper:
% "Group-Managed Real Options", by Lorenzo Garlappi, Ron Giammarino, and Ali Lazrak
% 
% June 29, 2021
%***********************************************************************************

% The figures in the appendices are reproduced and saved in the folder
% named 'Figures' under the current working directory

% The trading votes example is located at line 970

close all
clearvars

% Create folders 'Figures' and 'Results' if they do not exist
if ~exist('Figures', 'dir')
   mkdir('Figures')
end

if ~exist('Results', 'dir')
   mkdir('Results')
end

%% Which figures would you like to reproduce?

% Please put down 1 for figures that you
% would like to reproduce and 0 otherwise

% Supports any combinations of figures

all_figures = 1;
figure_1 = 0; 
figure_2 = 0;
figure_3 = 0;
figure_4 = 0;
figure_5 = 0;

% initialize values
[idx_figure_1, idx_figure_2, idx_figure_3, idx_figure_4, idx_figure_5] = deal([], [], [], [], []); 

% set get_figure indices (see line 97)
if all_figures == 1
    idx_figures = (1:13);
end
if figure_1 == 1
    idx_figure_1 = (1:3);
end
if figure_2 == 1
    idx_figure_2 = (13);
end
if figure_3 == 1
    idx_figure_3 = (4:11);
end
if figure_4 == 1
    idx_figure_4 = (4:11);
end
if figure_5 == 1
    idx_figure_5 = (12);
end

if all_figures == 0
    if sum(idx_figure_3) == sum(idx_figure_4)
        idx_figure_4 = [];
    end
    idx_figures = [idx_figure_1 idx_figure_3 idx_figure_4 idx_figure_5 idx_figure_2];
else 
    [figure_1, figure_2, figure_3, figure_4, figure_5] = deal(1, 1, 1, 1, 1); 
end

%% Parameters

r     = 0.05;
mu_OP = 4;
mu_P  = 0.01;
mu_O  = mu_OP * mu_P;
mu_M  = r/2;

mu_PP = 0;
mu_OO = mu_O + (r - mu_O)/2;

n_mu  = 100;
mu_M_vec  = linspace(mu_P, mu_O, n_mu);

scale = 1;

I = 1.0 * scale;

A_over_I_cases  = [0.05, 0.5, 0.6, 0.875, 0.9];  %Reversibility
sigma_cases     = [0.05, 0.3, 0.5];  %Volatility

LINEWIDTH = 2.0;
   
 %% Get Figures
 
for get_figure = idx_figures
    
    % What does 'get_figure' do? 
    % get_figure == 1: Figure 1 Panel A
    % get_figure == 2: Figure 1 Panel B
    % get_figure == 3: Figure 1 Panel C
    % get_figure == 4 to 9: Calculates results for Figure 3 & 4 (saved in the 'Results' folder)
    % get_figure == 10: Figure 3 Panel A & Figure 4 Panel A
    % get_figure == 11: Figure 3 Panel B & Figure 4 Panel B
    % get_figure == 12: Figure 5
    % get_figure == 13: Figure 2 (takes longer to run)
    
    if get_figure == 1
        idx_IA    = 3;
        idx_sigma = 2;
        PLOTS = 1;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 2
        idx_IA    = 4;
        idx_sigma = 2;
        PLOTS = 1;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 3
        idx_IA    = 3;
        idx_sigma = 3;
        PLOTS = 1;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 4
        idx_IA = 1;
        idx_sigma = 2;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 5
        idx_IA = 2;
        idx_sigma = 2;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 6
        idx_IA =5;
        idx_sigma = 2;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 7
        idx_IA = 1;
        idx_sigma = 3;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 8
        idx_IA = 2;
        idx_sigma = 3;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 9
        idx_IA = 5;
        idx_sigma = 3;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 10
        idx_IA = 1;
        idx_sigma = 1;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 1;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 11
        idx_IA = 1;
        idx_sigma = 1;
        PLOTS = 0;
        LOAD_RESULTS_sigma_high = 1;
        LOAD_RESULTS_sigma_med  = 0;
    elseif get_figure == 12
        idx_IA = 5;
        idx_sigma = 2;
        PLOTS = 1;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
    elseif get_figure == 13
        idx_IA = 1;
        idx_sigma = 1;
        PLOTS = 0;
        LOAD_RESULTS_sigma_med  = 0;
        LOAD_RESULTS_sigma_high = 0;
        
    end

    A     = A_over_I_cases(idx_IA) * I;
    sigma = sigma_cases(idx_sigma);
    
    if LOAD_RESULTS_sigma_med == 1
        sigma_val = 0.30;
        sigma_case = sprintf('sigma_med');
    elseif LOAD_RESULTS_sigma_high == 1
        sigma_val = 0.50;
        sigma_case = sprintf('sigma_high');
    end

    if LOAD_RESULTS_sigma_med ==1 || LOAD_RESULTS_sigma_high ==1
        result_file_IA_low  = sprintf('Results/results_A_over_I_0.05_sigma_%g.mat', sigma_val);
        result_file_IA_med  = sprintf('Results/results_A_over_I_0.5_sigma_%g.mat', sigma_val);
        result_file_IA_high = sprintf('Results/results_A_over_I_0.9_sigma_%g.mat', sigma_val);

        inertia_fig_file         = sprintf('Figures/Inertia_sigma_%g.eps', sigma_val);
        underinvestment_fig_file = sprintf('Figures/Underinvestment_sigma_%g.eps', sigma_val);
    end
    
    %% Individual Investment and Abandoment thresholds

    XA_PP = XA_ind_fcn(mu_PP, r, sigma, A);
    XA_P  = XA_ind_fcn(mu_P,  r, sigma, A);
    XA_M  = XA_ind_fcn(mu_M,  r, sigma, A);
    XA_O  = XA_ind_fcn(mu_O,  r, sigma, A);
    XA_OO = XA_ind_fcn(mu_OO, r, sigma, A);

    XA_M_vec = XA_ind_fcn(mu_M_vec, r, sigma, A);

    XI_PP = XI_ind_fcn(mu_PP, r, sigma, I, A, XA_PP);
    XI_P  = XI_ind_fcn(mu_P,  r, sigma, I, A, XA_P);
    XI_M  = XI_ind_fcn(mu_M,  r, sigma, I, A, XA_M);
    XI_O  = XI_ind_fcn(mu_O,  r, sigma, I, A, XA_O);
    XI_OO = XI_ind_fcn(mu_OO, r, sigma, I, A, XA_OO);


    XI_PP_SB_Maj_vec(1:n_mu) = NaN;
    XI_P_SB_Maj_vec(1:n_mu)  = NaN;
    XI_M_SB_Maj_vec(1:n_mu)  = NaN;
    XI_O_SB_Maj_vec(1:n_mu)  = NaN;
    XI_OO_SB_Maj_vec(1:n_mu) = NaN;


    XI_PP_SB_Una_vec(1:n_mu) = NaN;
    XI_P_SB_Una_vec(1:n_mu)  = NaN;
    XI_M_SB_Una_vec(1:n_mu)  = NaN;
    XI_O_SB_Una_vec(1:n_mu)  = NaN;
    XI_OO_SB_Una_vec(1:n_mu) = NaN;

    for i = 1:n_mu

        mu_M = mu_M_vec(i);

        mu_vec = [mu_PP, mu_P, mu_M, mu_O, mu_OO];

        XI_PP_SB_Maj_vec(i) = XI_group_fcn(mu_PP, mu_vec, r, sigma, I, A, 'Majority');
        XI_P_SB_Maj_vec(i)  = XI_group_fcn(mu_P,  mu_vec, r, sigma, I, A, 'Majority');
        XI_M_SB_Maj_vec(i)  = XI_group_fcn(mu_M,  mu_vec, r, sigma, I, A, 'Majority');
        XI_O_SB_Maj_vec(i)  = XI_group_fcn(mu_O,  mu_vec, r, sigma, I, A, 'Majority');
        XI_OO_SB_Maj_vec(i) = XI_group_fcn(mu_OO, mu_vec, r, sigma, I, A, 'Majority');

        XI_PP_SB_Una_vec(i) = XI_group_fcn(mu_PP, mu_vec, r, sigma, I, A, 'Unanimity');
        XI_P_SB_Una_vec(i)  = XI_group_fcn(mu_P,  mu_vec, r, sigma, I, A, 'Unanimity');
        XI_M_SB_Una_vec(i)  = XI_group_fcn(mu_M,  mu_vec, r, sigma, I, A, 'Unanimity');
        XI_O_SB_Una_vec(i)  = XI_group_fcn(mu_O,  mu_vec, r, sigma, I, A, 'Unanimity');
        XI_OO_SB_Una_vec(i) = XI_group_fcn(mu_OO, mu_vec, r, sigma, I, A, 'Unanimity');


    end

    %% Thresholds as a function of polarization (spread)


    mu_M_center = r/2;

    min_pol = 1e-4;
    max_pol = min(abs(0-mu_M_center), (r-min_pol - mu_M_center));

    n_pol = 1000;
    pol_vec = linspace(min_pol,max_pol,n_pol);


    XI_P_FB_pol(1:n_mu)  = NaN;
    XI_M_FB_pol(1:n_mu)  = NaN;
    XI_O_FB_pol(1:n_mu)  = NaN;

    XI_P_SB_Maj_pol(1:n_mu)  = NaN;
    XI_M_SB_Maj_pol(1:n_mu)  = NaN;
    XI_O_SB_Maj_pol(1:n_mu)  = NaN;

    XI_P_SB_Una_pol(1:n_mu)  = NaN;
    XI_M_SB_Una_pol(1:n_mu)  = NaN;
    XI_O_SB_Una_pol(1:n_mu)  = NaN;


    for i = 1:n_pol

        pol = pol_vec(i);
        mu_P_val = mu_M_center - pol;
        mu_M_val = mu_M_center;
        mu_O_val = mu_M_center + pol;

        mu_vec = [mu_P_val,mu_M_val,mu_O_val];

        % Individual policies
        XA_P_ind_pol = XA_ind_fcn(mu_P_val, r, sigma, A);
        XA_M_ind_pol = XA_ind_fcn(mu_M_val, r, sigma, A);
        XA_O_ind_pol = XA_ind_fcn(mu_O_val, r, sigma, A);

        XI_P_FB_pol(i)  = XI_ind_fcn(mu_P_val,  r, sigma, I, A, XA_P_ind_pol);
        XI_M_FB_pol(i)  = XI_ind_fcn(mu_M_val,  r, sigma, I, A, XA_M_ind_pol);
        XI_O_FB_pol(i)  = XI_ind_fcn(mu_O_val,  r, sigma, I, A, XA_O_ind_pol);



        % Group policies

        XI_P_SB_Maj_pol(i) = XI_group_fcn(mu_P_val,  mu_vec, r, sigma, I, A, 'Majority');
        XI_M_SB_Maj_pol(i) = XI_group_fcn(mu_M_val,  mu_vec, r, sigma, I, A, 'Majority');
        XI_O_SB_Maj_pol(i) = XI_group_fcn(mu_O_val,  mu_vec, r, sigma, I, A, 'Majority');

        XI_P_SB_Una_pol(i) = XI_group_fcn(mu_P_val,  mu_vec, r, sigma, I, A, 'Unanimity');
        XI_M_SB_Una_pol(i) = XI_group_fcn(mu_M_val,  mu_vec, r, sigma, I, A, 'Unanimity');
        XI_O_SB_Una_pol(i) = XI_group_fcn(mu_O_val,  mu_vec, r, sigma, I, A, 'Unanimity');

    end

    %% Inertia

    XI_FB     = [XI_P_FB_pol', XI_M_FB_pol',XI_O_FB_pol'];
    XI_SB_Maj = [XI_P_SB_Maj_pol', XI_M_SB_Maj_pol', XI_O_SB_Maj_pol'];
    XI_SB_Una = [XI_P_SB_Una_pol', XI_M_SB_Una_pol', XI_O_SB_Una_pol'];


    [~, n_Agents] = size(XI_SB_Maj);

    k_Maj = ceil(n_Agents/2);
    k_Una = n_Agents;

    k = k_Maj;

    [a_Maj,b_Maj] = mink([XI_P_SB_Maj_pol', XI_M_SB_Maj_pol',XI_O_SB_Maj_pol'],k_Maj,2);
    [a_Una,b_Una] = mink([XI_P_SB_Una_pol', XI_M_SB_Una_pol',XI_O_SB_Una_pol'],k_Una,2);

    XI_G_Majority  = a_Maj(:,k_Maj);
    idx_G_Majority = b_Maj(:,k_Maj);

    XI_G_Unanimity  = a_Una(:,k_Una);
    idx_G_Unanimity = b_Una(:,k_Una);

    Inertia_Majority  = XI_G_Majority ;
    Inertia_Unanimity = XI_G_Unanimity;


    XI_G_Maj_P_pivot(1:n_pol,1) = NaN;
    XI_G_Maj_M_pivot(1:n_pol,1) = NaN;
    XI_G_Maj_O_pivot(1:n_pol,1) = NaN;

    XI_G_Una_P_pivot(1:n_pol,1) = NaN;
    XI_G_Una_M_pivot(1:n_pol,1) = NaN;
    XI_G_Una_O_pivot(1:n_pol,1) = NaN;


    XI_G_Maj_P_pivot(idx_G_Majority==1) = XI_G_Majority(idx_G_Majority==1);
    XI_G_Maj_M_pivot(idx_G_Majority==2) = XI_G_Majority(idx_G_Majority==2);
    XI_G_Maj_O_pivot(idx_G_Majority==3) = XI_G_Majority(idx_G_Majority==3);

    XI_G_Una_P_pivot(idx_G_Unanimity==1) = XI_G_Unanimity(idx_G_Unanimity==1);
    XI_G_Una_M_pivot(idx_G_Unanimity==2) = XI_G_Unanimity(idx_G_Unanimity==2);
    XI_G_Una_O_pivot(idx_G_Unanimity==3) = XI_G_Unanimity(idx_G_Unanimity==3);


    %% Underinvestment--delta


    delta_Majority  = compute_delta(XI_G_Majority,  XI_M_FB_pol', XA_M_ind_pol, XA_M_ind_pol, mu_M_center, r, sigma, I, A);

    delta_Unanimity = NaN(n_pol,1);
    for i = 1:n_pol

        pol = pol_vec(i);
        mu_P_val = mu_M_center - pol;
        mu_M_val = mu_M_center;
        mu_O_val = mu_M_center + pol;

        XA_P_ind_pol = XA_ind_fcn(mu_P_val, r, sigma, A);
        XA_O_ind_pol = XA_ind_fcn(mu_O_val, r, sigma, A);


        delta_Unanimity(i) = compute_delta(XI_G_Unanimity(i), XI_P_FB_pol(i), XA_O_ind_pol, XA_P_ind_pol, mu_P_val, r, sigma, I, A);
    end



    delta_Maj_P_pivot(1:n_pol,1) = NaN;
    delta_Maj_M_pivot(1:n_pol,1) = NaN;
    delta_Maj_O_pivot(1:n_pol,1) = NaN;


    delta_Maj_P_pivot(idx_G_Majority==1) = delta_Majority(idx_G_Majority==1);
    delta_Maj_M_pivot(idx_G_Majority==2) = delta_Majority(idx_G_Majority==2);
    delta_Maj_O_pivot(idx_G_Majority==3) = delta_Majority(idx_G_Majority==3);



    delta_Una_P_pivot(1:n_pol,1) = NaN;
    delta_Una_M_pivot(1:n_pol,1) = NaN;
    delta_Una_O_pivot(1:n_pol,1) = NaN;

    delta_Una_P_pivot(idx_G_Unanimity==1) = delta_Unanimity(idx_G_Unanimity==1);
    delta_Una_M_pivot(idx_G_Unanimity==2) = delta_Unanimity(idx_G_Unanimity==2);
    delta_Una_O_pivot(idx_G_Unanimity==3) = delta_Unanimity(idx_G_Unanimity==3);


    %% Regions of parameters for which M is not pivotal

    if get_figure == 13
        
        mu_M_center = r/2;

        min_pol = 0;
        max_pol = min(abs(0-mu_M_center), (r-1e-5 - mu_M_center));
        n_pol   = 500;
        pol_cases = linspace(min_pol,max_pol,n_pol);

        min_AI   = 0;
        max_AI   = 1;
        n_A      = 500;
        AI_cases = linspace(min_AI, max_AI, n_A);
        A_cases = I * AI_cases;

        sigma_cases = [0.1:0.01:0.19, 0.2:0.1:1];
        n_sigma  = length(sigma_cases);

        XI_FB_Flag(1:n_pol, 1:n_A, 1:n_sigma)  = NaN;
        median_Flag(1:n_pol, 1:n_A, 1:n_sigma) = NaN;


        for i_pol = 1:n_pol
            disp(i_pol)
            pol   = pol_cases(i_pol);

            mu_P_val = mu_M_center - pol;
            mu_M_val = mu_M_center;
            mu_O_val = mu_M_center + pol;

            mu_vec = [mu_P_val,mu_M_val,mu_O_val];
            for i_A = 1:n_A

                A     = A_cases(i_A);

                for i_sigma = 1:n_sigma

                    sigma = sigma_cases(i_sigma);

                    XI_P_FB_pol  = XI_ind_fcn(mu_P_val,  r, sigma, I, A, XA_P_ind_pol);
                    XI_M_FB_pol  = XI_ind_fcn(mu_M_val,  r, sigma, I, A, XA_M_ind_pol);
                    XI_O_FB_pol  = XI_ind_fcn(mu_O_val,  r, sigma, I, A, XA_O_ind_pol);

                    if (XI_O_FB_pol < XI_M_FB_pol) && (XI_M_FB_pol < XI_P_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 0;

                    elseif (XI_O_FB_pol < XI_P_FB_pol) && (XI_P_FB_pol < XI_M_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 1;

                    elseif (XI_P_FB_pol < XI_M_FB_pol) && (XI_M_FB_pol < XI_O_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 2;

                    elseif (XI_P_FB_pol < XI_O_FB_pol) && (XI_O_FB_pol < XI_M_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 3;

                    elseif (XI_M_FB_pol < XI_O_FB_pol) && (XI_O_FB_pol < XI_P_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 4;

                    elseif (XI_M_FB_pol < XI_P_FB_pol) && (XI_P_FB_pol < XI_O_FB_pol)

                        XI_FB_Flag(i_pol, i_A, i_sigma) = 5;
                    end

                    % Group policies

                    XI_P_SB_Maj_pol = XI_group_fcn(mu_P_val,  mu_vec, r, sigma, I, A, 'Majority');
                    XI_M_SB_Maj_pol = XI_group_fcn(mu_M_val,  mu_vec, r, sigma, I, A, 'Majority');
                    XI_O_SB_Maj_pol = XI_group_fcn(mu_O_val,  mu_vec, r, sigma, I, A, 'Majority');

                    if (XI_M_SB_Maj_pol < XI_P_SB_Maj_pol && XI_P_SB_Maj_pol < XI_O_SB_Maj_pol) ||...
                            (XI_O_SB_Maj_pol < XI_P_SB_Maj_pol && XI_P_SB_Maj_pol < XI_M_SB_Maj_pol)

                        median_Flag(i_pol, i_A, i_sigma) = 1;  % P is pivotal

                    elseif (XI_M_SB_Maj_pol < XI_O_SB_Maj_pol && XI_O_SB_Maj_pol < XI_P_SB_Maj_pol) ||...
                            (XI_P_SB_Maj_pol < XI_O_SB_Maj_pol && XI_O_SB_Maj_pol < XI_M_SB_Maj_pol)

                        median_Flag(i_pol, i_A, i_sigma) = 2;  % O is pivotal

                    else

                        median_Flag(i_pol, i_A, i_sigma) = 0;  % M is pivotal
                    end



                end
            end
        end


        %% Plot heatmap
        close all
        mycolors = [.9 .9 .9; .8 .8 .9; .7 .7 .9];

        vol_case = 6;

        f = figure(1);
        colormap(f, mycolors)
        mesh(pol_cases, AI_cases,median_Flag(:,:,vol_case)', 'FaceColor','flat'); shg
        view([0,90])
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        ylabel('$A/I$','Interpreter','latex','fontsize',12)
        sigma_val = sigma_cases(vol_case);

        % Create textbox
        annotation(figure(1),'textbox',...
            [0.469642857142857 0.439476190476191 0.0848214285714285 0.0642857142857143],...
            'String',{'$n_I = M$'},'Interpreter','latex','fontsize',12);

        % Create textarrow
        annotation(figure(1),'textarrow',[0.626785714285714 0.833928571428571],...
            [0.829952380952381 0.892857142857143],'String',{'$n_I = O$'},'Interpreter','latex','fontsize',12);

        % Create textarrow
        annotation(figure(1),'textarrow',[0.742857142857143 0.896428571428571],...
        [0.646619047619048 0.890476190476191],'String',{'$n_I=P$'},'Interpreter','latex','fontsize',12);

         fig_file = sprintf('Figures/Region_AI_pol_sigma_%g.eps', sigma_val);
         print(fig_file, '-depsc2');    


         vol_case = 12;

        f = figure(2);
        colormap(f, mycolors)
        mesh(pol_cases, AI_cases,median_Flag(:,:,vol_case)', 'FaceColor','flat'); shg
        view([0,90])
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        ylabel('$A/I$','Interpreter','latex','fontsize',12)
        sigma_val = sigma_cases(vol_case);

        % Create textbox
        annotation(figure(2),'textbox',...
            [0.6875 0.77757142857143 0.0885791778564453 0.0523809523809524],...
            'String',{'$n_I =P$'},'Interpreter','latex','fontsize',12);

        % Create textbox
        annotation(figure(2),'textbox',...
            [0.585714285714284 0.651380952380956 0.0887515204293388 0.0523809523809524],...
            'String',{'$n_I = O$'},'Interpreter','latex','fontsize',12);

        % Create textbox
        annotation(figure(2),'textbox',...
            [0.389285714285713 0.370428571428577 0.0939052854265486 0.0523809523809524],...
            'String',{'$n_I = M$'},'Interpreter','latex','fontsize',12);

        fig_file = sprintf('Figures/Region_AI_pol_sigma_%g.eps', sigma_val);
         print(fig_file, '-depsc2');    



        vol_case = 14;

        f = figure(3);
        colormap(f, mycolors)
        mesh(pol_cases, AI_cases,median_Flag(:,:,vol_case)', 'FaceColor','flat'); shg
        view([0,90])
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        ylabel('$A/I$','Interpreter','latex','fontsize',12)
        sigma_val = sigma_cases(vol_case);


        % Create textbox
        annotation(figure(3),'textbox',...
            [0.716071428571428 0.637095238095239 0.0885791778564453 0.0523809523809524],...
            'String',{'$n_I = P $'},'Interpreter','latex','fontsize',12);

        % Create textbox
        annotation(figure(3),'textbox',...
            [0.530357142857143 0.491857142857144 0.0887515204293388 0.0523809523809524],...
            'String',{'$n_I = O$'},'Interpreter','latex','fontsize',12);

        % Create textbox
        annotation(figure(3),'textbox',...
            [0.364285714285713 0.251380952380955 0.0939052854265486 0.0523809523809524],...
            'String',{'$n_I = M$'},'Interpreter','latex','fontsize',12);

        fig_file = sprintf('Figures/Region_AI_pol_sigma_%g.eps', sigma_val);
         print(fig_file, '-depsc2');    

    end


    %% Save results
    
    SAVE_RESULTS = 0;
    if figure_3 == 1 && get_figure >= 4 && get_figure <= 9
        SAVE_RESULTS = 1;
    elseif figure_4 == 1 && get_figure >= 4 && get_figure <= 9
        SAVE_RESULTS = 1;
    elseif all_figures == 1 && get_figure >= 4 && get_figure <= 9
        SAVE_RESULTS = 1;
    end
    if SAVE_RESULTS == 1 
        
        outfile = sprintf('Results/results_A_over_I_%g_sigma_%g.mat', A_over_I_cases(idx_IA), sigma);
        save( outfile );

    end


    %% Plots

    if PLOTS ==1

        figure(1)
        plot(pol_vec, XI_P_FB_pol, 'b', pol_vec, XI_M_FB_pol, 'r', pol_vec, XI_O_FB_pol, 'g', 'linewidth', LINEWIDTH)
        legend('Pessimist, $X_P^{I,*}$','Median, $X_M^{I,*}$','Optimist, $X_O^{I,*}$', 'Location',  'Best')
        h = legend; set(h,'Interpreter','latex','fontsize',12); legend('boxoff');
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        shg

        if figure_1 == 1 && idx_IA == 3 && idx_sigma == 2
            fig_file_FB = sprintf('Figures/XI_FB_pol_A_over_I_%g_sigma_%g.eps', A_over_I_cases(idx_IA),sigma);
            print(fig_file_FB, '-depsc2')
        elseif figure_1 == 1 && idx_IA == 4 && idx_sigma == 2
            fig_file_FB = sprintf('Figures/XI_FB_pol_A_over_I_%g_sigma_%g.eps', A_over_I_cases(idx_IA),sigma);
            print(fig_file_FB, '-depsc2')
        elseif figure_1 == 1 && idx_IA == 3 && idx_sigma == 3
            fig_file_FB = sprintf('Figures/XI_FB_pol_A_over_I_%g_sigma_%g.eps', A_over_I_cases(idx_IA),sigma);
            print(fig_file_FB, '-depsc2')
        end

        figure(9)
        step = n_pol/100;
        plot(pol_vec(1:step:n_pol),XI_G_Maj_P_pivot(1:step:n_pol), 'b+',...
            pol_vec(1:step:n_pol), XI_G_Maj_M_pivot(1:step:n_pol), 'r+',...
            pol_vec(1:step/2:n_pol), XI_G_Maj_O_pivot(1:step/2:n_pol), 'g+',...
            pol_vec(1:step:n_pol), XI_G_Majority(1:step:n_pol), 'k:',...
            pol_vec(1:step:n_pol), XI_G_Una_P_pivot(1:step:n_pol), 'bo',...
            pol_vec(1:step:n_pol), XI_G_Una_M_pivot(1:step:n_pol), 'ro',...
            pol_vec(1:step/2:n_pol), XI_G_Una_O_pivot(1:step/2:n_pol), 'go',...
            pol_vec(1:step:n_pol), XI_G_Unanimity(1:step:n_pol), 'k:',...
            'linewidth', 1.2, 'MarkerSize',5)
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        ylabel('$X_G^I$','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])

        if idx_IA == 2
            annotation(figure(9),'textarrow',[0.414285714285714 0.516071428571428],...
                [0.639476190476191 0.504761904761905],'String',{'Unanimity'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(9),'textarrow',[0.696428571428571 0.783928571428571],...
                [0.418047619047619 0.342857142857143],'String',{'Majority'},'Interpreter','latex','fontsize',12);

        elseif idx_IA == 5
            annotation(figure(9),'textarrow',[0.517857142857143 0.594642857142857],...
                [0.510904761904762 0.380952380952381],'String',{'Unanimity'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(9),'textarrow',[0.716071428571429 0.635714285714286],...
                [0.270428571428571 0.152380952380952],'String',{'Majority'},'Interpreter','latex','fontsize',12);
        end        
        shg
        hold off

        if figure_5 == 1 && idx_IA == 5 && idx_sigma == 2
            fig_file = sprintf('Figures/XIG_pol_Maj_vs_Unanimity_A_over_I_%g_sigma_%g.eps', A_over_I_cases(idx_IA), sigma);
            print(fig_file, '-depsc2');
        end

        figure(10)
        step = n_pol/100;
        plot(pol_vec(1:step:n_pol),delta_Maj_P_pivot(1:step:n_pol)*100, 'b+',...
            pol_vec(1:step:n_pol),delta_Maj_M_pivot(1:step:n_pol)*100, 'r+',...
            pol_vec(1:step:n_pol),delta_Maj_O_pivot(1:step:n_pol)*100, 'g+',...
            pol_vec(1:step:n_pol),delta_Majority(1:step:n_pol)*100, 'k:',...
            pol_vec(1:step:n_pol),delta_Una_P_pivot(1:step:n_pol)*100, 'bo',...
            pol_vec(1:step:n_pol),delta_Una_M_pivot(1:step:n_pol)*100, 'ro',...
            pol_vec(1:step:n_pol),delta_Una_O_pivot(1:step:n_pol)*100, 'go',...
            pol_vec(1:step:n_pol),delta_Unanimity(1:step:n_pol)*100, 'k:',...
            'linewidth', 1.2, 'MarkerSize',5)
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        ylabel('$\delta$ (\%)','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        shg
        hold off

        if idx_IA == 2
            annotation(figure(10),'textarrow',[0.614285714285714 0.682142857142857],...
                [0.527571428571429 0.423809523809524],'String',{'Unanimity'},'Interpreter','latex','fontsize',12);
            % Create textarrow
            annotation(figure(10),'textarrow',[0.739285714285714 0.808928571428571],...
                [0.237095238095238 0.15],'String',{'Majority'},'Interpreter','latex','fontsize',12);
        elseif idx_IA == 5
            annotation(figure(10),'textarrow',[0.535714285714286 0.598214285714286],...
                [0.456142857142857 0.35],'String',{'Unanimity'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(10),'textarrow',[0.694642857142857 0.726785714285714],...
                [0.232333333333333 0.126190476190476],'String',{'Majority'},'Interpreter','latex','fontsize',12);

        end

        if figure_5 == 1 && idx_IA == 5 && idx_sigma == 2
            fig_file = sprintf('Figures/delta_pol_Maj_vs_Unanimity_A_over_I_%g_sigma_%g.eps', A_over_I_cases(idx_IA), sigma);
            print(fig_file, '-depsc2');
        end



    end
%%
    if LOAD_RESULTS_sigma_med ==1 || LOAD_RESULTS_sigma_high ==1
        

        load(result_file_IA_med)
        XI_G_Maj_P_pivot_AI_med  = XI_G_Maj_P_pivot;
        XI_G_Maj_M_pivot_AI_med  = XI_G_Maj_M_pivot;
        XI_G_Maj_O_pivot_AI_med  = XI_G_Maj_O_pivot;
        XI_G_Majority_AI_med     = XI_G_Majority;

        delta_Maj_P_pivot_AI_med = delta_Maj_P_pivot;
        delta_Maj_M_pivot_AI_med = delta_Maj_M_pivot;
        delta_Maj_O_pivot_AI_med = delta_Maj_O_pivot;
        delta_Majority_AI_med = delta_Majority;

        AI_med = A/I;

        load(result_file_IA_high)

        XI_G_Maj_P_pivot_AI_high  = XI_G_Maj_P_pivot;
        XI_G_Maj_M_pivot_AI_high  = XI_G_Maj_M_pivot;
        XI_G_Maj_O_pivot_AI_high  = XI_G_Maj_O_pivot;
        XI_G_Majority_AI_high     = XI_G_Majority;

        delta_Maj_P_pivot_AI_high = delta_Maj_P_pivot;
        delta_Maj_M_pivot_AI_high = delta_Maj_M_pivot;
        delta_Maj_O_pivot_AI_high = delta_Maj_O_pivot;
        delta_Majority_AI_high = delta_Majority;

        AI_high = A/I;


        load(result_file_IA_low)

        XI_G_Maj_P_pivot_AI_low  = XI_G_Maj_P_pivot;
        XI_G_Maj_M_pivot_AI_low  = XI_G_Maj_M_pivot;
        XI_G_Maj_O_pivot_AI_low  = XI_G_Maj_O_pivot;
        XI_G_Majority_AI_low     = XI_G_Majority;

        delta_Maj_P_pivot_AI_low = delta_Maj_P_pivot;
        delta_Maj_M_pivot_AI_low = delta_Maj_M_pivot;
        delta_Maj_O_pivot_AI_low = delta_Maj_O_pivot;
        delta_Majority_AI_low = delta_Majority;

        AI_low = A/I;



        step = n_pol/50;
        XI_mat = [XI_G_Maj_P_pivot_AI_med(1:step:n_pol),XI_G_Maj_M_pivot_AI_med(1:step:n_pol),XI_G_Maj_O_pivot_AI_med(1:step:n_pol),XI_G_Majority_AI_med(1:step:n_pol),...
            XI_G_Maj_P_pivot_AI_low(1:step:n_pol),XI_G_Maj_M_pivot_AI_low(1:step:n_pol),XI_G_Maj_O_pivot_AI_low(1:step:n_pol),XI_G_Majority_AI_low(1:step:n_pol),...
            XI_G_Maj_P_pivot_AI_high(1:step:n_pol),XI_G_Maj_M_pivot_AI_high(1:step:n_pol),XI_G_Maj_O_pivot_AI_high(1:step:n_pol),XI_G_Majority_AI_high(1:step:n_pol)];

        y_min = min(min(XI_mat));
        y_max = max(max(XI_mat));

        clf('reset')
        
        figure(1)
        plot(pol_vec(1:step:n_pol),XI_G_Maj_P_pivot_AI_low(1:step:n_pol), 'bo',...
            pol_vec(1:step:n_pol),XI_G_Maj_M_pivot_AI_low(1:step:n_pol), 'ro',...
            pol_vec(1:step:n_pol),XI_G_Maj_O_pivot_AI_low(1:step:n_pol), 'go',...
            pol_vec(1:step:n_pol),XI_G_Majority_AI_low(1:step:n_pol), 'k:',...
            pol_vec(1:step:n_pol),XI_G_Maj_P_pivot_AI_med(1:step:n_pol), 'b+',...
            pol_vec(1:step:n_pol),XI_G_Maj_M_pivot_AI_med(1:step:n_pol), 'r+',...
            pol_vec(1:step/5:n_pol),XI_G_Maj_O_pivot_AI_med(1:step/5:n_pol), 'g+',...
            pol_vec(1:step:n_pol),XI_G_Majority_AI_med(1:step:n_pol), 'k:',...
            pol_vec(1:step:n_pol),XI_G_Maj_P_pivot_AI_high(1:step:n_pol), 'b*',...
            pol_vec(1:step:n_pol),XI_G_Maj_M_pivot_AI_high(1:step:n_pol), 'r*',...
            pol_vec(1:step:n_pol),XI_G_Maj_O_pivot_AI_high(1:step:n_pol), 'g*',...
            pol_vec(1:step:n_pol),XI_G_Majority_AI_high(1:step:n_pol), 'k:',...
            'linewidth', 1.5, 'MarkerSize',5)
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        ylabel('$X_G^I$','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        ylim([y_min, y_max])

        if strcmp(sigma_case,'sigma_med') ==1

            % Create textarrow
            annotation(figure(1),'textarrow',[0.425 0.507142857142857],...
                [0.777571428571429 0.652380952380952],'String',{'Irreversible ($\circ$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(1),'textarrow',[0.269642857142857 0.432142857142857],...
                [0.246619047619048 0.142857142857143],'String',{'Reversible ($*$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(1),'textarrow',[0.646428571428571 0.539285714285714],...
                [0.349 0.457142857142857],'String',{'Baseline ($+$)'},'Interpreter','latex','fontsize',12);

        elseif strcmp(sigma_case,'sigma_high') ==1

            % Create textarrow
            annotation(figure(1),'textarrow',[0.353571428571429 0.425],...
                [0.80852380952381 0.911904761904762],'String',{'Irreversible ($\circ$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(1),'textarrow',[0.467857142857143 0.4],...
                [0.444238095238095 0.573809523809524],'String',{'Baseline ($+$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(1),'textarrow',[0.391071428571429 0.478571428571429],...
                [0.249 0.123809523809524],'String',{'Reversible ($*$)'},'Interpreter','latex','fontsize',12);

        end

        if figure_3 == 1
            inertia_fig_file = sprintf('Figures/Inertia_sigma_%g.eps', sigma);
            print(inertia_fig_file, '-depsc2');
        end


        clf('reset')
        
        figure(2)
        plot(pol_vec(1:step:n_pol),delta_Maj_P_pivot_AI_low(1:step:n_pol)*100, 'bo',...
            pol_vec(1:step:n_pol),delta_Maj_M_pivot_AI_low(1:step:n_pol)*100, 'ro',...
            pol_vec(1:step:n_pol),delta_Maj_O_pivot_AI_low(1:step:n_pol)*100, 'go',...
            pol_vec(1:step:n_pol),delta_Majority_AI_low(1:step:n_pol)*100, 'k:',...
            pol_vec(1:step/2:n_pol),delta_Maj_P_pivot_AI_med(1:step/2:n_pol)*100, 'b+',...
            pol_vec(1:step:n_pol),delta_Maj_M_pivot_AI_med(1:step:n_pol)*100, 'r+',...
            pol_vec(1:n_pol),delta_Maj_O_pivot_AI_med(1:n_pol)*100, 'g+',...
            pol_vec(1:step:n_pol),delta_Majority_AI_med(1:step:n_pol)*100, 'k:',...
            pol_vec(1:step:n_pol),delta_Maj_P_pivot_AI_high(1:step:n_pol)*100, 'b*',...
            pol_vec(1:step:n_pol),delta_Maj_M_pivot_AI_high(1:step:n_pol)*100, 'r*',...
            pol_vec(1:step:n_pol),delta_Maj_O_pivot_AI_high(1:step:n_pol)*100, 'g*',...
            pol_vec(1:step:n_pol),delta_Majority_AI_high(1:step:n_pol)*100, 'k:',...
            'linewidth', 1.5, 'MarkerSize',5)
        xlabel('Low~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~Polarization~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~High','Interpreter','latex','fontsize',12)
        ylabel('$\delta$ (\%)','Interpreter','latex','fontsize',12)
        set(gca,'xTick',[])
        shg

             
        if strcmp(sigma_case,'sigma_med') ==1

            % Create textarrow
            annotation(figure(2),'textarrow',[0.532142857142857 0.841071428571429],...
                [0.246619047619048 0.111904761904762],'String',{'Irreversible ($\circ$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(2),'textarrow',[0.667857142857143 0.880357142857143],...
                [0.327571428571429 0.169047619047619],'String',{'Reversible ($*$)'},'Interpreter','latex','fontsize',12);

            % Create textarrow
            annotation(figure(2),'textarrow',[0.6375 0.810714285714286],...
                [0.487095238095238 0.385714285714286],'String',{'Baseline ($+$)'},'Interpreter','latex','fontsize',12);


        elseif strcmp(sigma_case,'sigma_high') ==1

            % Create textarrow
            annotation(figure(2),'textarrow',[0.528571428571429 0.625],...
                [0.444238095238095 0.285714285714286],'String',{'Baseline ($+$)'},'Interpreter','latex','fontsize',12);

            % % Create textarrow
            annotation(figure(2),'textarrow',[0.3375 0.325],...
                [0.268047619047619 0.123809523809524],...
                'String',{'Reversible ($*$) \& Irrevertsible ($\circ$)'},'Interpreter','latex','fontsize',12);
        end
        
        if figure_4 == 1
            underinvestment_fig_file = sprintf('Figures/Underinvestment_sigma_%g.eps', sigma);
            print(underinvestment_fig_file, '-depsc2')
        end
     

    end
 
    
    
    close all

end


 %% Trading votes example

    example.mu_P    = 0.01;
    example.mu_M    = 0.02;
    example.mu_O    = 0.04;
    example.r       = 0.05;
    example.sigma   = 0.50;
    example.L       = 394.239;
    example.I       = 100;
    example.A       = 65;

    example.XA_P = XA_ind_fcn(example.mu_P, example.r, example.sigma, example.A);
    example.XA_M = XA_ind_fcn(example.mu_M, example.r, example.sigma, example.A);
    example.XA_O = XA_ind_fcn(example.mu_O, example.r, example.sigma, example.A);

    example.XI_P = XI_ind_fcn(example.mu_P, example.r, example.sigma, example.I, example.A, example.XA_P);
    example.XI_M = XI_ind_fcn(example.mu_M, example.r, example.sigma, example.I, example.A, example.XA_M);
    example.XI_O = XI_ind_fcn(example.mu_O, example.r, example.sigma, example.I, example.A, example.XA_O);

    example.XI_P_SB = XI_ind_fcn(example.mu_P, example.r, example.sigma, example.I, example.A, example.XA_M);
    example.XI_M_SB = XI_ind_fcn(example.mu_M, example.r, example.sigma, example.I, example.A, example.XA_M);
    example.XI_O_SB = XI_ind_fcn(example.mu_O, example.r, example.sigma, example.I, example.A, example.XA_M);

    example.XI_G = median([example.XI_P_SB, example.XI_M_SB,example.XI_O_SB]);
    example.XA_G = median([example.XA_P, example.XA_M,example.XA_O]);

    example.x_0     = 0.9 * example.XI_M;

    example.V_P_under_O = compute_V(example.x_0, example.XI_O, example.I, example.XA_O, example.A, example.mu_P, example.sigma, example.r);
    example.V_M_under_O = compute_V(example.x_0, example.XI_O, example.I, example.XA_O, example.A, example.mu_M, example.sigma, example.r);
    example.V_O_under_O = compute_V(example.x_0, example.XI_O, example.I, example.XA_O, example.A, example.mu_O, example.sigma, example.r);

    example.V_P_under_M = compute_V(example.x_0, example.XI_M, example.I, example.XA_M, example.A, example.mu_P, example.sigma, example.r);
    example.V_M_under_M = compute_V(example.x_0, example.XI_M, example.I, example.XA_M, example.A, example.mu_M, example.sigma, example.r);
    example.V_O_under_M = compute_V(example.x_0, example.XI_M, example.I, example.XA_M, example.A, example.mu_O, example.sigma, example.r);

    example.V_P_under_P = compute_V(example.x_0, example.XI_P, example.I, example.XA_P, example.A, example.mu_P, example.sigma, example.r);
    example.V_M_under_P = compute_V(example.x_0, example.XI_P, example.I, example.XA_P, example.A, example.mu_M, example.sigma, example.r);
    example.V_O_under_P = compute_V(example.x_0, example.XI_P, example.I, example.XA_P, example.A, example.mu_O, example.sigma, example.r);

    
    example.V_P_under_O_minus_L = example.V_P_under_O - example.L;
    example.V_M_under_O_minus_L = example.V_M_under_O - example.L;
    example.V_O_under_O_minus_L = example.V_O_under_O - example.L;

    example.V_P_under_M_minus_L = example.V_P_under_M - example.L;
    example.V_M_under_M_minus_L = example.V_M_under_M - example.L;
    example.V_O_under_M_minus_L = example.V_O_under_M - example.L;

    example.V_P_under_P_minus_L = example.V_P_under_P - example.L;
    example.V_M_under_P_minus_L = example.V_M_under_P - example.L;
    example.V_O_under_P_minus_L = example.V_O_under_P - example.L;
    


    example.V_P_solo  = compute_V(example.x_0, example.XI_P, example.I, example.XA_P, example.A, example.mu_P, example.sigma, example.r);
    example.V_M_solo  = compute_V(example.x_0, example.XI_M, example.I, example.XA_M, example.A, example.mu_M, example.sigma, example.r);
    example.V_O_solo  = compute_V(example.x_0, example.XI_O, example.I, example.XA_O, example.A, example.mu_O, example.sigma, example.r);
    example.V_M_group = compute_V(example.x_0, example.XI_G, example.I, example.XA_G, example.A, example.mu_M, example.sigma, example.r);

    example

%% Auxiliary Functions


function pi_val = compute_pi(x,X, mu, sigma,r)

m = compute_m(mu, sigma, r);
q = compute_q(mu, sigma, r);

pi_val =  (x/X).^m .* (x>=X) + (x/X).^q .* (x<X);
end

function W_val = compute_W(x, XA, A, mu, sigma, r)

pi_val = compute_pi(x,XA,mu, sigma, r);

W_val = (x./(r-mu) + (A-XA./(r-mu)) * pi_val) .* (x>=XA) + A * (x<XA);

end


function V_val = compute_V(x, XI, I, XA, A, mu, sigma, r)

W_x    = compute_W(x, XA, A, mu, sigma, r);
W_XI   = compute_W(XI, XA, A, mu, sigma, r);
pi_val = compute_pi(x, XI, mu, sigma, r);

V_val = (W_x - I).*(x>=XI) + (W_XI - I) .* pi_val .* (x<XI);


end


function XI_group_val = XI_group_fcn(mu, mu_vec, r, sigma, I, A, rule)
% This function determines the second-best investment threshold for group
% of any odd-size and for Majority and/or Unanimity

n_Members  = length(mu_vec);

if strcmp(rule, 'Majority')
    mu_Pivot_A = mu_vec(ceil(n_Members/2));
elseif strcmp(rule, 'Unanimity')
    mu_Pivot_A = max(mu_vec);
end

XA = XA_ind_fcn(mu_Pivot_A, r, sigma, A);


XI_group_val = XI_ind_fcn(mu, r, sigma, I, A, XA);

end

function XI_ind_val = XI_ind_fcn(mu, r, sigma, I, A, XA)

m = compute_m(mu, sigma, r);
q = compute_q(mu, sigma, r);


XI_guess =q/(q-1) * (r-mu) * I;

[XI_ind_val, fval, exitflag] = fzero(@(x) FOC_fcn(x, mu, r, m, q, I, A, XA), XI_guess);

% This is to make sure that solutions are found for mu very close to r
if exitflag ~=1
    MAXITER = 10;
    iter = 0;
    while exitflag ~=1 && iter < MAXITER
        
        XI_guess = XI_guess + .1;
        iter = iter +1;
        [XI_ind_val, fval, exitflag] = fzero(@(x) FOC_fcn(x, mu, r, m, q, I, A, XA), XI_guess);
        
    end
end

end



function FOC_val = FOC_fcn(x, mu, r, m, q, I, A, XA)

FOC_val = x ./ (r-mu)*(q - 1) + (A - XA/(r-mu)) .* (x/XA).^m * (q - m) - q * I;


end


function XA_ind_val = XA_ind_fcn(mu, r, sigma, A)

m = compute_m(mu, sigma, r);


XA_ind_val = m./(m-1) * A .* (r - mu);

end


function m_val = compute_m(mu, sigma, r)

mu_minus_pt5_sig2 = mu - 0.5*sigma^2;

m_val = 1/sigma^2* ( -(mu_minus_pt5_sig2) - ( mu_minus_pt5_sig2.^2 + 2*sigma^2*r ).^(0.5) );

end

function q_val = compute_q(mu, sigma, r)

mu_minus_pt5_sig2 = mu - 0.5*sigma^2;

q_val = 1/sigma^2* ( -(mu_minus_pt5_sig2) + ( mu_minus_pt5_sig2.^2 + 2*sigma^2*r ).^(0.5) );

end



function delta_val = compute_delta(XI_Group, XI_Solo, XA_Group, XA_Solo, mu, r, sigma, I, A)

Gamma_Group_val = compute_W(XI_Group, XA_Group, A, mu, sigma, r) - I;
Gamma_Solo_val  = compute_W(XI_Solo,  XA_Solo, A, mu, sigma, r)  - I;

q = compute_q(mu, sigma, r);

delta_val = (Gamma_Solo_val./Gamma_Group_val).^(1./q) .* (XI_Group./XI_Solo) - 1;


end